% This code reports variations in rate stability cost estimates by
% aggregate risks. 
function analyze_rstability(cost_rs, cost2_1, M, K, D_unmerge)

% total cost 
major=D_unmerge.major; % NNx 1
pr_s = D_unmerge.pr_s;

% - all states including zero increase states
Ecost_rs = sum(pr_s.*cost_rs,2);
Ecost_rs_pe = Ecost_rs./D_unmerge.s1; % per-enrollee cost

% - excluding zero increase states
pr_s_inc = pr_s(:, [2:5 7:10 12:15 17:20]);
cum_pr = sum(pr_s_inc,2);
cum_pr = repmat(cum_pr, 1, 16);
pr_s_inc = pr_s_inc./cum_pr;

cost_rs_inc = cost_rs(:, [2:5 7:10 12:15 17:20]);
Ecost_rs_inc = sum(pr_s_inc.*cost_rs_inc,2);
Ecost_rs_inc_pe = Ecost_rs_inc./D_unmerge.s1; % per-enrollee cost

% mean cost2_1 conditional on mm
Ecost2_1_mm = zeros(size(cost2_1,1),M);
for mm=1:M
    col1 = 1+K*(mm-1);
    col2 = K*mm;
    Ecost2_1_mm(:,mm) = sum(cost2_1(:,col1:col2).*D_unmerge.pr_k,2);
end

% mean cost2_1 conditional on k
Ecost2_1_k = zeros(size(cost2_1,1),K);
for kk=1:K
    Ecost2_1_k(:,kk) = sum(cost2_1(:,kk:K:end).*D_unmerge.pr_mm,2);
end

% for each interest rate value, integrate over NH shock
% - deno is zero if the observation never experices high interest rate in t=2
deno_rhigh = D_unmerge.pr_mm(:,1) + D_unmerge.pr_mm(:,2);
pr_temp = D_unmerge.pr_mm(:,1)./deno_rhigh;
Ecost2_1_rhigh_all = Ecost2_1_mm(:,1).*pr_temp + Ecost2_1_mm(:,2).*(1-pr_temp);
Ecost2_1_rhigh     = Ecost2_1_rhigh_all(deno_rhigh>0);

deno_rlow = D_unmerge.pr_mm(:,3) + D_unmerge.pr_mm(:,4);
pr_temp = D_unmerge.pr_mm(:,3)./deno_rlow;
Ecost2_1_rlow_all = Ecost2_1_mm(:,3).*pr_temp + Ecost2_1_mm(:,4).*(1-pr_temp);
Ecost2_1_rlow     = Ecost2_1_rlow_all(deno_rlow>0);

% for each NH value, integrate over interest shock
deno_nhlow = D_unmerge.pr_mm(:,1) + D_unmerge.pr_mm(:,3);
pr_temp = D_unmerge.pr_mm(:,1)./deno_nhlow;
Ecost2_1_nhlow_all = Ecost2_1_mm(:,1).*pr_temp + Ecost2_1_mm(:,3).*(1-pr_temp);
Ecost2_1_nhlow     = Ecost2_1_nhlow_all(deno_nhlow>0);

deno_nhhigh = D_unmerge.pr_mm(:,2) + D_unmerge.pr_mm(:,4);
pr_temp = D_unmerge.pr_mm(:,2)./deno_nhhigh;
Ecost2_1_nhhigh_all = Ecost2_1_mm(:,2).*pr_temp + Ecost2_1_mm(:,4).*(1-pr_temp);
Ecost2_1_nhhigh     = Ecost2_1_nhhigh_all(deno_nhhigh>0);

Ecost2_1_kbottom = 0.5*Ecost2_1_k(:,2) + 0.5*Ecost2_1_k(:,3);
Ecost2_1_ktop = 0.5*Ecost2_1_k(:,4) + 0.5*Ecost2_1_k(:,5);

disp(' ')
disp('Mean cost2_1 ')
disp(' ')
disp(': By the residual agg risk')
disp(['- conditional on being in bottom half = ', num2str(mean(Ecost2_1_kbottom))])
disp(['- conditional on being in top half = ', num2str(mean(Ecost2_1_ktop))])

disp(' ')
disp(': By NH cost shock')
disp(['- conditional on NH cost low       = ', num2str(mean(Ecost2_1_nhlow))])
disp(['- conditional on when NH cost high = ', num2str(mean(Ecost2_1_nhhigh))])

disp(' ')
disp(': By interest rate shock')
disp(['- conditional on interest rate high = ', num2str(mean(Ecost2_1_rhigh))])
disp(['- conditional on interest rate low  = ', num2str(mean(Ecost2_1_rlow))])
